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Abstract 



We consider divergence form elliptic operators in dimension n > 2 with 
L°° coefficients. Although solutions of these operators are only Holder con- 
tinuous, we show that they are differentiable (C 1,Q ) with respect to harmonic 
coordinates. It follows that numerical homogenization can be extended to 
situations where the medium has no ergodicity at small scales and is char- 
acterized by a continuum of scales by transferring a new metric in addition 
to traditional averaged (homogenized) quantities from subgrid scales into 
computational scales and error bounds can be given. This numerical ho- 
mogenization method can also be used as a compression tool for differential 
operators. 

1 Introduction and main results 

Let Q be a bounded and convex domain of class C 2 . We consider the 
following benchmark PDE 



Where g is a function in L°°(Q). And x — > a(x) is a mapping from Q, 
into the space of positive definite symmetric matrices. We assume a to be 
symmetric, uniformly elliptic with entries in L°°(0,). 
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• Is it possible to up-scale 

Homogenization theory (|15j. |52j ) allows us to do so by transferring bulk 
(averaged) information from sub- grid scales into computational scales. This 
transfer from a numerical homogenization point of view is justified under 
two fundamental assumptions: 

• Ergodicity at small scales, and 

• scale separation. 

Can we get rid of these assumptions? 

• Can we numerically homogenize (jl.lj) when a is arbitrary? in particu- 
lar when a is characterized by a continuum of scales with no ergodicity 
at small scales? 

What do we mean by numerical homogenization when a is arbitrary? It is 
important to recall that F. Murat and L. Tartar's theory of H-convergence 
jtilj provides a mathematical framework for analysis of composites in com- 
plete generality, without any need for geometrical hypotheses such as pe- 
riodicity or randomness. This theory is based on a powerful tool called 
compensated compactness or div-curl lemma introduced in the 70's by Mu- 
rat and Tartar [HO], [ZHj which has been characterized by a wide range of 
applications and refinements |25| . Here we consider homogenization from a 
slightly different point of view: we want to solve (jl.lj) on a coarse mesh and 
we want to understand which information should be transferred from fine 
scales to coarse scales when the entries of a are arbitrary. For that purpose 
we need a new form of compensation given in section 11.11 

It is important to observe that if one needs to solve (jl.lj) only one time 
(with one g) the method proposed here does not reduce the number of nu- 
merical operations 1 . Indeed we need to compute (jl.lj) n-times (n being the 
space dimension) with in the right hand side in (jl.lj) and linear boundary 
conditions. However if one needs to solve (jl.lj) for a large number (M) of 
different right hand sides g (M » n, which would happen if one tries to 
optimize specific properties of u with respect to g) then the methods pro- 
posed here have a practical use since they basically say that after solving 
(jl.lj) n times it is sufficient to solve (jl.lj) on a coarse mesh (with TV 0,01 nodes 
instead of N for instance) M times. 

Let us recall that fast methods based on hierarchical matrices are avail- 
able [13 EES D3B Ej for solving JDJ in 0(N (In N) n+3 ) operations. 

Finally the point of view of this paper is to observe the "redundancy" 
of solutions of (jl.lj) at small scales rather than "fast computation" . Indeed 
it is not obvious that (jl.lj) can be homogenized when the medium does 

1 It will be shown in that for parabolic and reaction-diffusion equations the situation 
is different: the methodology introduced in this paper can be used to reduce the number of 
operations even when one needs to solve these equations only once. 
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not satisfy the usual periodicity, ergodicity or scale-separation assumptions. 
Moreover if (jl.ljl can indeed be homogenized it is important to understand 
what minimal quantity of information should be kept from small scales to 
obtain an accurate homogenized operator? Once the correct coarse parame- 
ters are identified (the bulk quantities and the up-scaled metric) , one can try 
to model, estimate or simulate them but the first step is to identify them. 



1.1 A new form of compensation 

To introduce the new form of compensation, we need to introduce the so 
called a-harmonic coordinates associated to i.e. the weak solution of 

the following boundary value problem 



div aVF = in tt 
Fix) = x on <9f2. 



11.2) 



By p.2|) we mean that F is a n-dimensional vector field F(x) = (iq(x), . . . , F n (x)) 
such that each of its entries satisfies 



div aS/Fi = in tt 
Fi(x) = Xi on dQ. 



(1.3) 



The new compensation phenomenon is controlled by the following object: 
Definition 1.1. We write 

a := *VFaVF. (1.4) 
We write \i a the anisotropic distortion of a defined by 

/Vax(CT(x))x 

H„ := esssupseo I > ' J ■ (1.5) 

Where A max (M) (A m j n (M)) denote the maximal (minimal) eigenvalue of M. 

Definition 1.2. In dimension n = 2, we say that a is stable if and only if 
fjL a < oo and there exist a constant e > such that (Trace(er)) e G L 1 (Q) 

Remark 1.1. According to [H in dimension two if a is smooth then a is 
stable. According to [JQ, F is always an homeomorphism in dimension two 
even with a^j 6 L°°(Q). 

Theorem 1.1. Assume that a is stable and n = 2. Then there exist con- 
stants a > and C > such that (VF) _1 Vu G C Q (0) and 

|| (VF)- 1 ^^^ < C||<7|| Loo(n) . (1.6) 
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Remark 1.2. The constant a depends on 0, A max (a)/A m i n (a) and fj, a . The 
constant C depends on the constants above, A m i n (a) and || ( Trace(cr)) llii(n)" 
We use the notation A max (a) := sup xgf2 supi£i =1 *£a£ and 
Amax(a) := inf^eQ mf| €[=1 *£a£. 

Remark 1.3. If one considers a sequence a e such that /v e and 
||(Trace(a e )) _ ~ e || ri/ -, are uniformly bounded away from oo and A m i n (a e ) 
and A max (a e ) are uniformly bounded away from and oo then Ijl.Gj) is uni- 
formly true. If we consider a periodically oscillating sequence a e (x) = a(f ), 
/Vi < oo and || ( Trace (o"i)) 1 ^H^i^) < 00 then (|1.6|) is uniformly true. 

Remark 1.4. It is easy to check from the proof of ll.ll that if ( Trace(<7 e )) 1 G 
L q (£l) (with g > 2) then it is possible to replace by ||<7||z,p(n) in 

(|1.6j) with p depending on q. More precisely if (Trace(<7 e )) 1 € L°°(Q) then 
it is possible to replace ||<7||z,°o(n) by 1 1 | j z, 2 += (o) m Q1-6JI for any e > 0. 
Remark 1.5. We don't need a to be symmetric to obtain fTTTl but for the 
clarity of the presentation we have chosen to restrict ourselves to that case. 
We have not analyzed Neumann boundary conditions, this will done in a 
further work. 

Remark 1.6. We will use the notation V_pn := (Vi ? )~ 1 Vti. In dimension 
two, it is known ([I], jS], pQ) that the determinant of V-F is strictly positive 
almost everywhere and the object V^u is well defined. In dimension three 
and higher V fu is well defined when a is stable. 



This phenomenon can be observed numerically. In figures 1(a) and 2 (a 
a is given by a product of random functions oscillating over a continuum of 
scales. The entries of the matrix VF are in L p (figure |2(b)| , the entries of 
the gradient of u in the Euclidean metric are in L p (figures |l(b) and 2(c) ) 
yet (VF) _1 Vu is Holder continuous (figures [T(c)1 and [2(d)] ). 



Let us now introduce the compensation phenomenon in dimension n > 3. 
We call Pa- the Cordes parameter associated to a defined by 

/ ( Trace[cr(x)]) 2 \ 
Pa := es SS up xeQ [n - ^———j . (1.7) 



Observe that since P a is also given by 

a ( ( Yli=l \,<t(x)) \ ,.. „x 

Pa = esssup xgn (n = — -= J. (1.8) 

where (\% % m ) denotes the eigenvalues of M, it is a measure of the anisotropy 
of a. 

Definition 1.3. In dimension n > 3, we say that a is stable if and only 
if, Pa- < 1 and if n < 4 that there exist a constant e > such that 
(Trace(cj)) 2 G L l (VL) 
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(a) a in log scale. (b) one of the entries VF. 




(c) Vit. (d) (VFj-'Vu. 

Figure 2: Change of metric on the torus. 
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Remark 1.7. According to ^] and [2j| in dimension three and higher a can 
be unstable even if a is smooth. We refer to figureHOfor an explicit example. 

Let us write 

an i 
(j2 my?)* d X y (i.9) 
- ! t,i=i 



Theorem 1.2. Assume that a is stable and n > 3. T/ien F is an auto- 
morphism on Q. Moreover there exist constants p > 2 and C > suc/i i/ioi 
uoF- 1 e W^ p (fl) and 

\\ uoF ~ 1 \\w£' p (n) - C \\9\\l™(q,)- (1-10) 

Remark 1.8. The constant » depends on n,0 and /3o-. The constant C de- 
pends on the constants above, A m j n (a) and if n < 4 on || ( Trace(cr)) 2 1 1 x, 1 (O) " 

In the following theorem we do not need to assume to be convex. 

Theorem 1.3. Assume n > 2 and (Trace(a)) -1 G L°°(Q). Let p > 2. 
There exist a constant C* = C*(n,d£l) > such that if (3 a < C* then there 
exists a real number 7 > depending only on n,Q and p such that 

IKVFj^VtiH^dt^CIHIi,,^. (1.11) 

Remark 1.9. The constant C in (|1.11|) depends on n, 7, 0, C*, A m i n (a), 
\\a\\L°°(n), and ||(Trace(cr)) _1 || ioo(nT) . 



1.2 Dimensionality reduction 

Observe that (|1.1|) is a priori an infinite dimensional problem since a and g 
can be irregular at all scales. Yet according to theorem ll.ll and ll.21 whatever 
the choice of g, at small scales, solutions to (f 1 . If) are correlated to F which 
lives in a functional space of dimension n. More precisely we will propose a 
rigorous justification of a variation of the multi-scale finite element method 2 
introduced by Hou and Wu 49 in its form refined by Allaire and Brizzi 
in situations where the medium is not assumed to be periodic or ergodic 
(these methods are already rigorously justified when the medium is periodic 

021.0). 

Let 7/j be a coarse conformal mesh on Q composed of n-simplices (triangles 
in dimension two and tetrahedra in dimension three). Here h is the usual 
resolution of the mesh defined as the maximal length of the edges of the 
tessellation. Let us call 7(7^) the maximum over of the n-simplices K of 7^ 

2 Let us recall that the multi-scale finite element method is inspired from Tartar's oscillating 
test functions [23] 
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of the ratio between the radius of the smallest ball containing K and the 
largest ball inscribed in K. We assume "f{Th) to be uniformly bounded in 
h. 

We write Vh C H l (yL) the set of piecewise linear functions on the coarse 
mesh vanishing at the boundary of the tessellation. We write Mh the set 
of interior nodes of the tessellation and <pi (i G Mh) the usual nodal basis 
function of Vh satisfying 

i Pi{yj) = ^ij- (1-12) 

We consider the elements (ipi)iej\f h defined by 



ipi := (fii o F(x). 



(1.13) 





(a) ipi (b) ipi = ipi o F 

Figure 3: The Galerkin elements 



Let us write Uh the solution of the Galerkin scheme associated to (|1.2|) 
based on the elements {ipi)i^f h - Observe that the number of elements is on 
the order of h~ n and we have the following theorem 

Theorem 1.4. Assume that a is stable and n = 2. Then there exist con- 
stants a, C > such that 



UhWii 1 < C/i a ||g|| LO o(Q). 



(1.14) 



Remark 1.10. The constant a depends only on n,Cl and [i a . The constant 
C depends on the objects mentioned above plus A m i n (a), A max (a), "f{Th) and 
||(Trace(a))~ 1-e || il(n) . 

Remark 1.11. Theorem 11.41 is also valid with a = 1 as in theorem 11.51 The 
only difference between these two theorems lies in the constant C. In the 
proof theorem 11.41 we use the property u o F^ 1 € C ,Q! (0) and in the proof 
of theorem ll.5l we use the property u o F^ 1 £ W 2,2 (£l). 
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Remark 1.12. Let us recall that Uh is an element of the space Xh spanned 
by (ipi)iej\f h obtained as the solution of the linear problem 



Where a[, ] is the bilinear form on Hq(Q) defined by 



(1.15) 



a[v,w] := I 
Jn 



VvaVw. 



(1.16) 



It follows from theorem 11.41 that solutions to IJl.lj) live in the iT -norm 
neighborhood of a low dimensional space. 

Remark 1.13. The proof of theorem ll.4l is done for the exact function tpi and 
not for its discrete version on a fine mesh. If a is regular at a given small 
scale ho then it is easy to check that theorem 11.41 remains valid as long as the 
edges of the fine mesh are smaller than Iiq. A more intriguing case is when 
a is discrete and discontinuous on a fine mesh. Numerical experiments show 
that theorems such as 1 1 . 41 and 1 1 . 1 1 remain valid but to justify them for the 
discrete version of the harmonic coordinates F and elements one would 
have to adapt our theorems to the discrete setting. In order to remain 
concise we have chosen to not include that adaptation in this paper. 

Remark 1.14. We keep the composition rule used in The only difference 
between the elements Ijl.lMJI and the ones proposed by Hou, Wu, Allaire and 
Brizzi lies in the fact that we use the global solution to (jl.2|) and not a local 
one computed on each triangle of the coarse mesh through an over-sampling 
technique. We refer to remark T3. II for further comments. 

Remark 1.15. Write S the stiffness matrix a[ipi,ipj]. S^ 1 is in general dense 
(characterized by N 2 entries where N is the number of nodes of the mesh) . 
Yet surprisingly by combining theorem 11.41 with theorem 5.4 of ^3] one 
can obtain that S^ 1 can be approximated (in L 2 norm) with hierarchi- 
cal matrix Mu such that the matrix vector product by Mh requires only 
0(N(inN) n+3 ) operations. 

In dimension n > 3 we have the following estimate. 
Theorem 1.5. Assume that a is stable, n > 3 and 

|| ( Trace(cr)) £,00(01 ^ 00 ' Then there exist constants a,C > such that 



Remark 1.16. The constant C depends on n, 7(^/1), ^2, A max (a), A m j n (a) 
and ||(Trace(cr)) _1 || ioo(or 



u - u h \\m < Ch\\g\\ L00 (ay 



(1.17) 
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-0.1 ■ 
-0.2 ■ 

-0.6 -0.4 -0.2 0.2 

Figure 4: Support of the elements ipi and fy. 



1.3 Galerkin with localized elements 

For the clarity of the paper we will restrict ourselves from now on to di- 
mension two, the generalization of the statements to higher dimensions is 
conditioned on the stability of a (and the application of theorem II, 3j) . 

The elements (jl.l3|) can be highly distorted and non local (figure 0J) since 

support(^j) := F~ l (support^)) . (1-18) 

It follows that the elements ipi are piecewise linear on a fine mesh different 
from the one on which a is defined and F has been computed. Is it possible 
to avoid that difficulty by solving on a coarse mesh with localized 

elements? The answer is yes but the price to pay for the localization will 
be the discontinuity of the elements and the fact that the accuracy of the 
method will depend on a weak aspect ratio of the triangles of the tessellation 
in the metric induced by F. More precisely consider a triangle K of the 
tessellation, call a, b, c the nodes of K and 9 the interior angle of the triangle 
(F(a), F(b), F(c)) which is the closest to tt/2. We call weak aspect ratio of 
the triangle K in the metric induced by F the quantity 

»»<»< A '> : = (1 ' 19> 

So r]^i n (K) is large if the triangle (F(a), F(b), F(c)) is flat (all its interior 
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angles are close to or n). We define 



Vmin ■= SUp 1]mm( K )- (l- 2 0) 

Ker h 

Let us recall that although the coefficients of the PDE are irregular it 
is well known 7Q a that F is Holder continuous. Thus it makes sense to look 
at the value of F at a specific point. Now let v be a function defined on the 
nodes of the triangle K € Th, let us write a, b, c the nodes of that triangle. 
It is usual to look at the coarse gradient of v evaluated at the nodes of the 
triangle K, i.e. the vector defined by 

b — a \ 1 / v(b) — v(a) 



Vv(K) :=( . (1.21) 

\ c — a J \ v(c) — v(a) J 

If i]^ in (K) < oo then the following object called the gradient of v evaluated 
on the coarse mesh with respect to the metric induced by F is well defined. 

V P » W :=(^-£iin?|-f>). (1-22) 
V F(c) - F(a) J \ v{c) - v(a) J 

Definition 1.4. We say that the tessellation Th is not unadapted to F if 
and only if the determinant of VF(K) is strictly positive for all K 6 7^. 

Remark 1.17. Observe that if the tessellation 7^ is not unadapted to F then 
r]^ in (K) < oo, the definition II .41 contains the additional condition that there 
is no inversion in the images of the triangles of Th by F. 

Now consider the nodal elements {S,i)i^M h i defined by 

= k (L23) 

I Vi?^(x) = constant within each K € Th. 

If the mesh is not unadapted to F then the elements (figure EJ) (|1.23|) are 
well defined and given by 

J^.( x ) = l+ (F(x) -F(xi))V F tpi(K) if i~K and xeK 
1 &(0) = i n other cases. 

(1-24) 

Where the notation i ~ K means that i is a node of K. Observe that 
the elements £j are discontinuous at the boundaries of the triangles of the 
coarse mesh however they are easier to implement since they are localized 
in these triangles. Write Zh the vector space spanned by the functions £j. 
For K £ Th we write ax the bilinear form on H l (K) defined by 

a K [v,w] := / *VwzVw. (1.25) 
Jk 
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Figure 5: Localized Galerkin elements £j 



We will write i? 1 (7/ l ) the space of functions v 6 L 2 (ft) such that the re- 
striction of v to each triangle K belongs to H l (K). We will write for 
v t w£H\T h ) 

a*[v,w] := ^2 a K [v,w]. (1.26) 

K&T h 

The localized finite element method can be formulated in the following way: 
look for u* € Zh such that for all i E Afh, 

a*fo,uf] = &,g) LHQ) . (1.27) 

We have the following estimate 

Theorem 1.6. Assume that a is stable and that the mesh is not unadapted 
to F. Then there exists a constant a > such that 

{a*[u-uf\)* < Crj* min h a \\g\\ Loo{n) . (1-28) 

Remark 1.18. For a bilinear from B[,] we write B[v] := B[v,v]. 

Remark 1.19. The constant a depends only on n, fi, e and \i a . The constant 
C depends on the objects mentioned above plus || (Trace(er)J Hxifnv 
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Remark 1.20. The bilinear operator a*[, ] on Zh is characterized by a con- 
stant matrix within each triangle K £ 7/j equal to 

*(VF(iO) -1 ('VFaVF}^(VF(iO) (1-29) 

where < v >k means the average of v over K with respect to the Lebesgue 
measure (< v >k'-= voi(.ft:) Ik v ( x } d%-> Vol (if) being the Lebesgue measure 
of volume of K). 

The error bound given in theorem 11.61 is given in the norm induced by 
a*[.]. We would like to obtain an error bound with respect to the usual H 1 
norm. Observe that is discontinuous at the boundaries of the triangles of 
the coarse mesh so we have to find an accurate way to interpolate u$ in the 
whole space using its values at the nodes of the coarse mesh. Let us write 
F{J\fh) the image of the nodes of 7^ by F. Let us write T F the triangulation 
of F(Nh)- Let ipf be the standard piecewise linear nodal basis of T F . Let 
us write Jh the interpolation operator from the space of functions defined 
on the nodes of 7^ into H l {Q) defined by 

J h v{x) := Y, v ( x i)Vi oF ( x )- ( L3 °) 

Observe that for i £ Afh, v{xi) = Jhvixi)- We have the following estimate 

Theorem 1.7. Assume that a is stable and that the mesh is not unadapted 
to F. Then there exist constants a,Cf > such that 

\\u - Jhu f \\m(n) < Cfh a \\g\\L°°(n)- (1-31) 

Remark 1.21. The constant a depends only on n, Q and \i a . The constant 
Cj can be written 

C f := C^ in (min (C ax r^ ax , u*)j 2 (1.32) 

where C depends on the objects mentioned above plus A m i n (a), A max (a) and 
|| ( Trace(cr)) 1 e || L i(Q)- ??max is defined by where 9 the interior angle 
of the triangles of closest to or tt. ^ is defined by where 7 the 
interior angle of the triangles of T F closest to or tt. Moreover 

Vol(K F ) 

where K F is the triangle whose nodes are the images of the nodes of K by 
F. 
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1.4 Numerical homogenization from the informa- 
tion point of view. 

The Galerkin scheme described in 11.21 and 11.31 are based on elements con- 
taining the whole fine scale structure of F. This represents too much in- 
formation. We can wonder what minimal quantity of information should 
be kept from the scales in order to up-scale Q1.1JI ? We would like to keep 
an accurate version of (|1.1|) with minimal computer memory. This point 
touches the compression issue. Images can be compressed. Can the same 
thing be done with operators? 

This question has received an answer within the context of the fast mul- 
tiplication of vectors with fully populated special matrices arising in var- 
ious applications (35] • Let us recall the fast multipole method and the 
hierarchical multipole method designed by L. Greengard and V. Rokhlin 
|43| . Wavelet based methods have been designed by G. Beylkin, R. Coif- 
man and V. Rokhlin [HI I17[ I16j . The concept of Hierarchical matrices has 
been developed by W. Hackbusch et al. .43,. More precisely we refer to 
|13l |H1 1121 11UI II 1 j . The Hierarchical matrix method is based on a compres- 
sion of the inverse of the stiffness matrix (see remark ll.l5j) . Here we consider 
compression from the point of view of numerical homogenization. We look 
at the operator (|1.1|) as a bilinear form on Hq(Q) and we will use Vh as 
space of test functions to zoom at the operator associated to a at a given 
arbitrary resolution. 



The up-scaled or compressed operator, written Una will naturally be a bilin- 
ear form on the space of piecewise linear functions on the coarse mesh with 
Dirichlet boundary condition. 



The question is how to choose W^a? To answer that question we can integrate 
Ql-ljl against a test function in Vh, then we obtain that 



We will use the test function </> to "look at" the operator 1)1. lj) at the given 
resolution h. We can decompose the first term in the integral above as a 
sum of integrals over the triangles of the coarse mesh to obtain (we assume 
that a is stable) 




(1.34) 




(1.35) 




(1.36) 




(1.37) 
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Now V</> is constant within each triangle K € T^. (VF(x)^ 1 \7u(x) is 
Holder continuous thus almost a constant within each triangle K and equal 
to the gradient of u evaluated on the coarse mesh with respect to the metric 
induced by F, i.e. the following vector 



F(c) - F(a) J \ u(c) - u(a) 

Where a, b, c are the nodes of the triangle K. It follows that the tensor 
aVF can be averaged over each triangle of the coarse mesh and we will 
write < aVF >k its average. In conclusion a good candidate for the up- 
scaled operator U^a is the bilinear form given by the following formula: for 

v,w eVh 

U h a[v,w] := V f t Vv(aVF) K (VF(K)y 1 Ww. (1.39) 

K€T h jK 

Observe that the only information kept from the small scales in the com- 
pressed operator (jl.39j) are the bulk quantities (aVF)^ and the non aver- 
aged quantities F(b) — F(a) where a and b are nodes of the triangles of the 
coarse mesh. The latter quantity can be interpreted as a deformation of the 
coarse mesh induced by the small scales (or a new distance defining coarse 
gradients). In the particular case where a = M(^) and M is ergodic then 
as e | < aVF >k converges to the usual effective conductivity obtained 
from homogenization theory and X7F(K) converges to the identity matrix. 
It follows that the object (|1.39|) recovers the formulae obtained from ho- 
mogenization theory when the medium is ergodic and characterized by scale 
separation. Let us now show that formula p. 39)1 can be accurate beyond 
these assumptions. 

To estimate the compression accuracy we have to use the up-scaled op- 
erator lA^a to obtain an approximation of the linear interpolation of u on 
the coarse mesh. We look for u m £ Vh such that for all i € Mh, 

U h a[pi,u m ] = (<Pi,g) L 2( n) . (1.40) 

The price to pay for the loss of information on the small scales is the loss of 
ellipticity. This loss can be caused by two correlated factors: 

• The new metric can generate flat triangles. 

• The up-scaled operator can become singular. 

The first factor is due to the localization of the scheme. The second factor 
does not appear with Galerkin schemes. It is not observed in dimension 
two but it can't be avoided in dimension higher or equal to three in the 
sense that the up-scaled operator has no reason to remain elliptic and local. 
Indeed consider a box of dimension three, and set in that box empty tubes 



14 



of low boundary conductivity as shown in figure ||3 Set the left side of the 
box to temperature 0°c and the right side to temperature 100°c. Then an 
inversion in the temperature profile is produced around the critical points 
shown in figure El (see [I] and > instead of increasing from left to right 
in these region temperature decreases). Now as the operator is up-scaled, 
the information on the geometry of the tubes is lost but the inversion phe- 
nomenon remains in the loss of ellipticity and locality of the operator. We 
will address this issue further in a forthcoming paper. 



Critical points 




Figure 6: a in dimension three 



Nevertheless it is possible to prove that once stability is achieved then 
the method is accurate (if a is stable). More precisely, for a nodal function 
v let us define the homogeneous Dirichlet form on the graph induced by T^: 

ZJl«i-t>il 2 (1-41) 

We write i ~ j when those nodes share an edge on the coarse mesh. Let us 
define the following stability parameter of the scheme 

cm ■ f U h a[v,w] 

S := mi sup ^ r . (1-42) 

">ev»vev h (£ h [ v ])2(E h [ w ])2 

Observe that S m depends only on the up-scaled parameters so we have a 
control on the stability. 

Definition 1.5. We say that the scheme is stable if and only if S m > 0. 

Remark 1.22. Let us recall that for v G Vh, £h[ v ] can be bounded from below 
and above by the L 2 -norm of the gradient of v. More precisely 

£h[v] < ||V«||! 2 (n) < r] max £ h [v]. (1.43) 



r) max = 1/ sin(6 l ) where 6 is the closest interior angle of the triangles of Th 
to or 7T. 
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Remark 1.23. In practice in dimension two the condition number of the 
scheme associated to the up-scaled operator is as good as the one obtained 
from a Galerkin scheme by solving a local cell problem. 

Let us write X^u the linear interpolation of u over 7^: 

l h u := u(xi)(pi(x). (1.44) 

We have the following estimate 

Theorem 1.8. Assume that a and the scheme are stable and that the mesh 
is not unadapted to F. Then there exist constants ct,C m > such that 

\\l h u - u m \\ H i (n) < C m h a \\g\\ Lx{n) . (1.45) 

Remark 1.24. The constant a depends only on n,Q and [i a . The constant 
C m can be written 

C m ■= C^p. (1.46) 

where C depends on the objects mentioned above plus A m i n (a), A max (a) and 
||( Traced)) _1_e || Ll 




(a) u (b) u m (c) Triangle of the 

coarse mesh. 

Figure 7: u estimated with the up-scaled operator. 

The compressed operator allows us to capture the solution of (jl.l|) on a 
coarse mesh (figure[7J) . What information should be added to the compressed 
operator in order to obtain fine resolution approximation of u! The answer 
is a finer resolution of F (figure EJ). Indeed let be the interpolation 
operator introduced in (jl.HOj) . we then have the following estimate 
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Tiylor Expansion af Lie 




(a) Refined F. . (b) Refined F. (c) J h u. 

Figure 8: Taylor expansion with respect to the new metric. 



Theorem 1.9. Assume that a and the scheme are stable and that the mesh 
is not unadapted to F. Then there exist constants a, C m > such that 

h - JhU m \\ H i m < C m h a \\g\\ L oo {n) . (1-47) 

Remark 1.25. The constant a depends only on n,0 and \i a . The constant 
C m can be written 

C m :=c( V ™f™ x f (1.48) 

where C depends on the objects mentioned above plus A m j n (a), A max (a) and 
|| ( Trace(cr))~ 1 ~ e || il(0) . 

1.4.1 Coherent multi-resolution. 

We want to compress a physical system from a fine scale description (F) to 
a coarse scale description (C) there are two ways of doing so: 

• Either we up-scale directly from (F) to (C) 

• Either we do so in two steps: from (F) to an intermediate scale (I), 
then from (I) to (C). 

Now if the scales (F), (I) and (C) are not completely separated a technique 
solely based on averaging (up-scaling of bulk quantities) would produce two 
different results (depending on the presence of an intermediate step or not). 
Thus it is important to check the consistency of the numerical homogeniza- 
tion method if the metric information F(xi) — F(xj) is up-scaled in addition 
to traditional bulk quantities < aVF >k- 

Let T° , . . . ,T n be a multi-resolution tessellation of 0. Each T % is a 
regular conformal tessellation of £1. Moreover 7~ l+1 is a refinement of T l . 
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Let us write V 1 the space of piecewise linear functions on T l . We write B % 
the space of bilinear operators on T l . We want to compress (up-scale) the 
bilinear operator a[, ] on the multi-grid T°,...,T ra . We assume that the 
smallest scale n is fine enough to capture the irregularities of a, in that case 
we define a n such that for all v, u G V n 

a n [v,u] := a[v,u]. (1.49) 

Since the gradient of an element of V n is constant within each triangle of 
T n , a n [, ] can be defined by a mapping from T n onto M. n the space of n x n 
constant matrices. We will write a n (K) the constant matrix associated 3 to 
K G T n . Similarly each bilinear operator of B % can be defined by mapping 
from T l onto A4 n . We define for k < p, U k ' p the up-scaling operator mapping 
B p onto B k in the following way (we assume that the tessellations % are not 
unadapted to F and that the respective schemes are stable). Let B G B p . 

• Let F G V p be the solution of 

(B\v,F]=0 for all v G V p , 

P r* (L50) 

\r=x on 1 . 

where we have written T l the boundary of T\ 

• The bilinear form U k,p B is defined by its matrices lA k,p B(K) for K G 
T k . 

U k ' p B{K) := (BVF) K (VF(K))~ 1 . (1.51) 
(.)„ stands for the averaging operator 

^ VF ^ := Voi7M S Vol(T)B(T)VF(T). (1.52) 
For k < p < q, U satisfies the semi-group property 

Note also that U q,q = Id- In particular if we define for k G {1, . . . , n} 

a k := U k ' n a n (1.54) 

as the up-scaled operator, the following semi-group property is satisfied: for 

k < p < n 

a k := U k ' p a p . (1.55) 

Semi-group properties (|1.53f) . ()1.55j) are essential to the consistency and 
coherence of the numerical homogenization method. 



3 the bilinear form (|1.49|) can be written as a sum of integrals over K G T n . Vv and Vu are 
constant over these triangles, thus bilinear form over V n is determined by a matrix. 
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1.5 Numerical homogenization from the transport 
point of view. 

The elliptic operator appearing in (jl.lj) can be seen as the generator of a 
stochastic differential equation. This stochastic differential equation can re- 
flect the transport process of a pollutant in a highly heterogeneous medium 
such as soil. The following operator A — VW whose numerical homogeniza- 
tion is similar to the one of can represent a physical system evolving in 
a highly irregular energy landscape V. The simple fact that this evolution 
taking place in a continuous domain can be captured by a Markov chain 
evolving on a graph is far from being obvious |72| . Our point of view here 
is to accurately simulate a Markov chain living on a fine graph by an "up- 
scaled" Markov chain living on a coarse graph. The main question is how 
to choose the the jump rate 7^ of the random walk between the nodes of 
the coarse graph? The answer to that question is given by a finite volume 
method. 

Let us write 7^* the dual mesh associated to % t . can be obtained by 
drawing segments from the midpoints of the edges of the triangles of 7^ to 
an interior point in these triangles (the circumcenter to obtain a Voronoi 
tessellation but one can also choose the barycenter). 

Let us write Vi the control volume associated to the node i of the primal 
mesh and %i t ne characteristic function of Vi. The finite volume method 
can be expressed in the following way: look for u v £ Zh (Zh being the space 
spanned by the elements & introduced in Q1.23JI ) such that for all i € Mh 



Again, it follows from equation 1)1. 56|) that the only information kept from 
the scales are the usual bulk quantities (effective conductivities at the edges 
of the dual mesh) plus the metric information F(b) — F(a) where a and b 
are nodes of the triangles of the primal mesh. Observe also that it is possi- 
ble to generate with this finite volume method a coherent multi-resolution 
compression similar to the one introduced in subsection II .4. 11 According to 
()1.56j) the good choice for the jump rates of the random walk should be 



To properly describe the transport process one should look at a parabolic 
operator instead of the elliptic one. This issue will be addressed in |65| . 
we will restrict ourselves to the elliptic case characterizing the equilibrium 
properties of the random walk. Let us write S v the stability of the up-scaled 
finite volume operator. It is defined by 




(1.56) 



Hi = a *[Xi,£j] if i ~ 3 and % / j. 



(1.57) 



S 



w£Zhvey h (£ h [v})2(£ h [w})2 




(1.58) 
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Theorem 1.10. Assume that a and the scheme are stable (S v > 0) and 
that the mesh is not unadapted to F. Then there exist constants a,C v > 
such that 



and 



\\l h u - l h u v \\ H i {n) < C V) ih a \\g\\ Loa ^ (1.59) 



\u - J h u v \\ H T.{Q) < C Vt 2h a \\g\\ L oo^y (1.60) 



Remark 1.26. The constant a depends only on n,Cl and [i a . The constant 
C v i can be written 

C v ,i:=C^0^{X max (a))K (1.61) 
The constant C v % can be written 

C depends on the objects mentioned above plus A m i n (a), A max (a) and 
||( Traced)) _1_e || il(n) . 

Remark 1.27. Observe that we need the additional condition A max (c) < 00 
to prove the convergence of the method. Numerical experiments show that 
although the finite volume method keeps very few information from small 
scales it is more stable and accurate than the method presented in 11.41 (it is 
also more stable and almost as accurate as Galerkin method in which the 
whole fine scale scale structure of F is up-scaled). That is why we believe 
that the constants (|1.61|) and (|1.62|) are not optimal. 

1.6 Explicit formulae in laminar cases. 

The harmonic coordinates can be explicitly computed in dimension one. In 
this subsection we will analyze a toy model to understand the effect of the 
new metric when the coefficients of the partial differential equations are 
characterized by an infinite number of overlapping scales. Our point is to 
show that new metric becomes multi-fractal. 

Let Q := (0,1). Let V G L°°(J7) and write w(V) the weak solution of 
the following Dirichlet problem: 

f-Ie^div(e-^ V ^=/ 
I w = on dVL, 
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With / G L°°(0). We write w(V) the solution of Write /i + and //_ 

the probability measures defined on the Borelian subset of (0, 1) by 

u+0,x:=^ and u_0,d:=^ . (1.64) 

Define D(V) := (jje 2 ^ e'^^dz)' 1 . D(V) is to put into corre- 
spondence with the bulk quantity < aV-F > which is a constant in dimension 
one. Observe that (TTHal) can be explicitly rewritten 

In that sense the metric based numerical homogenization is exact in dimen- 
sion one. Equations such as (|1.65|) have been studied in dimension one in 
and [AH\ in order to introduce a measure-theoretic way of defining dif- 
ferential operators on fractal sets on the real line. 

Let V n be a sequence in L°°(0). We call a probability measure non degen- 
erate if and only if it is atom-less and the mass of any non void open subset 
of Q is strictly positive. We have the following theorem 

Theorem 1.11. If //T n and pYr weakly converge to non degenerate probabil- 
ity measures p + and p- then D(y n )w(V n ) converges pointwise to the unique 
solution of the following differential equation with Dirichlet boundary con- 
dition 

^ ^ ^ tp f (1 66) 

2 dp,- dp+ 

In the case of classical homogenization observe that p + and are simple 
Lebesgue probability measures on [0,1]. 

Write T the torus of dimension one and side of length one and let U G 
C X (T). Let p E N/{0, 1}. Write T p the scaling operator denied on the space 
of functions by T p U{x) := U{px). Write 

n-l 

S n U :=J2t p pU. (1.67) 

p=0 

Take V n = S n U in theorem II .111 Then by Perron- Frobenius-Ruelle theorem 
|76| . pY^ and weakly converge to some probability measure //+ and //_ 
(eigenvectors of the Ruelle transfer operator) and it is easy to check that 
they are non degenerate. Let us notice that similarly it is possible to show 
that — ^lnD(V n ) converges to the sum of topological pressures of U and 
—U with respect to the shift induced by the multiplication by p on the 
space of p-adic decompositions |66| . |14j . Theorem ll.lll is telling us that the 
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regularity of ip corresponds to the regularity of /i+[0, x] thus it is natural to 
wonder what is the regularity of that harmonic measure. To answer that 
question we will consider the paradigm of binomial measures. We refer to 
|HH] for a detailed introduction to this subject, for the sake of completeness 
we will recall its main lines below in our framework. We take U(x) G L°°(T) 
with (a ^ b) 

(a for < x < i 
2 U(x) = { k ~ 2 (1.68) 

1 b for ~ < x < 1. 



Let us write 



e b 



ran = r and mi = r. (1.69) 

e a + e b e a + e b 

Then /xT n weakly converges to //+ and that measure is self-similar in the 
sense that it satisfies 

fi+([x,y}) = m fi + ([2x,2y})+m lf i + ([2x -l,2y- 1]). (1.70) 

For x € (0, 1) we write x±X2 ■ ■ ■ its 2-adic decomposition in the sense that 
x = Y^=o x n2 n ■ We write I yit ... )yp the cylinder of x G (0,1) such that 
xi . . . x p = y\ ... y p . fi + is atom-less but singular with respect to the 
Lebesgue measure. Moreover it is easy to check [B5] that fi + and there- 
fore ip are Holder continuous and their Holder continuity exponent is not a 
constant. More precisely let us write 

rami 

and 

a(x) := lim a n (x). (1-72) 

n—»oo 

Whenever this limit exists. This limit exists for almost all x (with respect 
to Lebesgue measure) and its value depends on the dyadic expansion of x. 
Writing l n (x) the number of one appearing in the first n digits of x we have 

a(x) = lim -(1 - ^M) log 2 (m ) - ^ Iog 2 (mi). (1.73) 
rwoo inn mn 

Thus a{x) can take all the values between — log 2 (mo) and — log 2 (mi). How- 
ever for almost all x with respect to Lebesgue measure 

a(x) = ~log 2 (m omi ) = ~{^- ~ He a + e b )). (1.74) 

Now it is possible to obtain from large deviation theory (j^J theorem 2 and 
ES|) that 

a n (x) E (a — e,a + e)j —t 1 + c*(a) (1-75) 
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P being the uniform probability measure on (0,1), c(q) = 1 — log 2 (?rtQ + 
ml) and c* denotes the Legendre transform of c, i.e. c*(a) = inf g (qa — 
c(q)) . Thus the metric associated to our up-scaling method is multi-fractal 
|4()j . |57j . Let us recall that multi- fractal formalism was originally 

introduced to describe the regularity of the velocity field in Turbulence jlU] 
and explaining intermittency. 



1.7 Literature. 

The issue of numerical homogenization partial differential equations with 
heterogeneous coefficients has received a great deal of attention and many 
methods have been proposed 4 . Let us mention a few of them. 

• Multi-scale finite element methods [3U], [M], EH], H3, EH, 02], El- 

• Multi-scale finite volume methods |54| . 

• Heterogeneous Multi-scale Methods 29 . 

• Wavelet based homogenization g3],|2B|,|231,D3|.|Z|.D21- 

• Residual free bubbles methods [2*0] . 

• Discontinuous enrichment methods [HI], |33| . 

• Partition of Unity Methods [2Z] • 

• Energy Minimizing Multi-grid Methods |75| . 

The methods mentioned above are part of a larger quest aimed at captur- 
ing high dimensional problems with a few coarse parameters [HH], EH [62| . 
|44j . Paraphrasing the outcome of a recent DOE workshop .21], we may 
understand the physics of multi-scale structures at each individual scale nev- 
ertheless without the capability to bridge the scales, a significant number of 
important scientific and engineering problems will remain out of reach. 



2 Proofs. 

2.1 Compensation. 

Let us prove theorem ll.il We need a variation of Campanato's result |22j on 
non-divergence form elliptic operators. Let us write for a symmetric matrix 
M, 

7i ,_ J2i=i Km (r> ,s 
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Justified in dimension one, in the case of periodic or ergodic media with scale separation or 



in the case of partial differential equations with sufficiently smooth coefficients 
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We consider the following Dirichlet problem: 

L M v = f (2.2) 
with Lm '■= Ya,j=i Mij(x)d{dj. We assume M to be elliptic and symmetric. 

Theorem 2.1. Assume that (3m < 1- If is convex, then, there exists a 
real number p > 2 depending only on n,Q and Pu such that if f £ L P (Q) 
the Dirichlet problem 12.2]) has a unique solution satisfying 



C 

\ W ^P(Q) < —\\vMf\\LP{Q)- (2-3) 

1 - R 2 
1 Pm 



Remark 2.1. Pm is the Cordes parameter (|1.8[) associated to M. 

Proof. Theorem |23 is a straightforward adaptation of theorem 1.2.1 of 56 , 
for the sake of completeness we will give the main lines of ideas leading to 
estimate Q2.3|l . Let us recall the Miranda- Talenti estimate |56| . 

Lemma 2.1. Let Q C R n be a bounded and convex domain of class C 2 . 
Then for each v € Wq ,2 (Q) it results 

[ V (didjv) 2 dx < [(Avfdx. (2.4) 
Jn ij=1 Jn 

The Laplacian A : Wq' p {Q) -> U>(fi) is an isomorphism for each p > 1. 
Let A^ 1 (p) be the inverse operator A -1 : L P (Q) — > Wq' p . It is clear from 
(|2.4j) that |[A~ 1 (2)|[ < 1. Let r E (2,oo), by the convexity of the norms we 
have 

IIA^CpJH <C(p) (2.5) 

with 

C{p) :=\\A~ l {r)\\^) . (2.6) 

Let v be a solution of (|2.2[l (we refer to |56j for the existence of v which is 
obtained from a fix point theorem), we have 

\M w ^ { n) ^ \\^ 1 {p)\\\\^ v \\lv(q)- (2.7) 

Observing that Av = vuf + Av — vmLmv one can obtain 

\\Av\\ LP{n) < \\vMf\\up(n) + \\^ v ~ vmLmvWlp^). (2.8) 

Then following the proof of theorem theorem 1.2.1 of jEB] we have 

r. n 

\\Av - u M L M v\\ P LP{n) < / tf( Y, dx - ( 2 - 9 ) 

Jn i,j=i 
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Let us choose p > 2 such that l-C(p)/?^ 2 > (l-/3^ 2 )/2. Combining (|2~?jl . 
(|2~%j) and (fZHj) we obtain that 

/ ( dx )" ^ ^%\Hd\\ m ny (2-10) 

Jn i,j=i 1 ~Pm 

Which leads to estimate (j2.3j) , □ 
Remembering the Sobolev embedding inequality 

H^Ha 1 -^) ^ c W v Ww*»m (2 - n) 

theorem 12. II implies the Holder continuity of v in dimension n = 2. 

We assume that a is stable. We write F~ l the inverse of F (which is 
well defined if a is stable). Let us write Q the symmetric positive matrix 
given by the following equation 

W : = ( , wnJ ° F ~ (y)- ( 2 - 12 ) 



I det (VF , 

Let us write u; the solution of the following equation: for all <p € Cg°, 
n ^ 

V QndidjW = j—^- -. (2.13) 

.4^ J J det(VF) oF-i v ; 

Let us now prove the following theorem, 

Theorem 2.2. Assume that a is stable and that £1 is convex. Then there 
exists constants p > 2, C > such that the solution of \2.13\) belongs to 
Wq' p (Q) and satisfies 

C 

ll v Hlw/ 2 ' p (c) < rlMU~(n)- (2-14) 

1-/33 



Remark 2.2. a depends on f2 and (3 a . C depends on A m i n (a) and if n < 4 
on || (Trace(a)) 7_2_ l L i (n) - 

Proof. Now let us observe that 

_ Trace(a) 



det(VF) o F- 1 Trace(o- 2 

Using the change of variables y = F(x) and choosing 1/q' + 1/q = 1 we 
obtain that 

ii U Q n ii Trace(o") . M 

" d^F ^ - ll 5 ll^'wllTr^p) (det(VF))P9 ll^^)- 

(2.16) 
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It is easy to check that 
ll^™*,^ < / n („)i- (2,7) 

For 2 < n < 4 we choose q = 1 in ()2.17j) for n > 5 we choose q = Then 
a direct application of theorem 12.11 and estimate (|2.15j) to equation Q2.1HJI 
implies the theorem (observe that Pq = f3 a ). □ 

Let (p € Cq°(Q). Write (p := tpo F^ 1 . Using theorem 12 .21 we obtain that 



<p, 



Qndidjw) =-(0,-, t—^ t) ■ (2-18) 

' J J T.2(a\ V' Uot UP nP-l/r.2m) v ' 



£ 2 (fi) V' |det(VF)| of-Vi'tfi) 



Using the change of variable y = F{x) we deduce that 

n 

^a ijmw) oF) =-^, 9 ) (2.19) 
Let us observe that 

n 

(Tij(didjw) o F = div (aVF((Vw) o F) V (2.20) 

i,i=l 

It follows after an integration by parts that (and observing that VF(Vw) o 
F = V(woF)) 

4ft l " oJf ] = (ft9)i J (fi) ( 2 - 21 ) 

It follows from the uniqueness of the solution of the Dirichlet problem (|2.21f) 
that w o F = u. Theorem 11.21 is then a straightforward consequence of 
theorem 12.21 and the equality u o F^ 1 = w. 

Remark 2.3. Using the change of variables y = F(x) we obtain for all ip G 

Co°°(^)) 

a[tp,u]=Q[tp,u] (2.22) 

The comparison between (|2.18|) and (|2.22j) indicates that Q is a divergence 
free matrix. We have not used that property of Q explicitly in our proof 



above but it is present implicitly in the deduction of (|2.21|) from (|2.20|) . 

Remark 2.4. The only place where we use the convexity of f2 is for the 
validity of lemma 12.11 (we refer to jSH] ) ■ 

The following lemma is a well known result obtained from the De Giorgi- 
Moser-Nash theory (j2Hl) |59j . ) of divergence form elliptic operators with 
discontinuous coefficients (more precisely we refer to |70| for the Global 
Holder regularity) 
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Lemma 2.2. There exists C, a' > depending on f2 and A max (a)/A m i n (a) 
such that F is a' Holder continuous and 

\\F\\ Ca/ < C. (2.23) 

Theorem 1 1.1 1 is a straightforward consequence of the Sobolev embedding 
inequality (|2.11|) . theorem 12.21 lemma l2~T1 and the fact that V^u = VuoF. 
Let us observe that in dimension two, we have 

1 1, 1 X , 

(/V + — ). (2.24) 



1 - p ff 2 ^ 
And the condition (3 a < 1 is equivalent to \i a < oo. 

2.1.1 Holder continuity for n > 3 or non-convexity of Q. 

In this subsection we will not assume Q, to be convex. Let N p ' (p,) (1 < p < 
oo, < A < n) be the weighted Morrey space formed by functions v : O — > R 
such that ^([^^(f^) < oo with 

IMIiVP.*(n) = SU P ( / \x - x \~ x \v(x)\ p ) P . (2.25) 

To obtain the Holder continuity of u o i* 1-1 in dimension n > 3 we will use 
corollary 4.1 of We will give the result of S. Leonardi below in a form 
adapted to our context. Consider the Dirichlet problem (|2.2|) . We do not 
assume to be bounded. We write W 2 ' P ' X (Q) the functions in W 2 ' P (Q) such 
that their second order derivatives are in NP'^^lj). 

Theorem 2.3. There exist a constant C* = C*(n,p, \,d£l) > such that 
if Pm < C* and f S N P ' X (Q) then the Dirichlet problem has a unique 

solution in W 2 ' P ' X n Wq ,p (Q). Moreover, ifO<X<n<p then S/v G C a (Q) 
with a = 1 — n/p and 

C 

l|Vu||o«(n) < x , . M J I/IUyA(n) ( 2 - 26 ) 

where C = C(n,p, A, 90). 

Theorem 11.31 is a straightforward application of theorem 12.31 

2.2 Dimensionality reduction. 

Let us prove theorem 11.41 We write Xh the linear space spanned by the 
elements ipi. The solution of the Galerkin scheme satisfies a[u — u^, v] = 
for all v € Xh. Thus 

a[u — Uh] = inf a[u — u^, u — v\. (2.27) 

vex h 
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It follows by Cauchy-Schwartz inequality that 



a[u — Ufi) < inf a[u — v}. 

vex h 



(2.28) 



Now writing v := v o F 1 and using the change of variable y = F(x) we 
obtain 

a[u - v] = Q[u - v}. (2.29) 
It follows in dimension n = 2 that 

D 



\u - u h \\ H i < 



Amin(a) w€V h 



inf II Vu — Vu>H ? 



L°°(n) 



with 



D := Trace 



/ 



VFaVF 



(2.30) 



(2.31) 



Thus using the following standard approximation properties of the elements 
ipi (see for instance |32]). 



inf ||Vu — Vt^H^/Q) < C^y('77 i )/i CK ||€t H^i, 



wev h 



cg' a (n) 



(2.32) 



we obtain that 



||« - < 7Cft)( > ^ J ^||Vu|| ca fe a . (2.33) 

We conclude by observing that for I € M n 

/¥VFaVFZ = inf / ' *(| + V f)a(l + V '/). (2.34) 

Theorem 1 1 . 41 b ecomes a direct consequence of (|2.3U[) and theorem 12.21 
In dimension n > 3 we obtain from (12.291) that 



i- ^ -^max(Q) . t — * _ 1,2 , .-, or , 

\\u - u h \\ m < j-r- mf Vu - Vfl) U n) . (2.35) 

Amin(o) ™eV A ; 

It is easy to obtain that 

A m ax(Q) < (det(a))^|(Trace(a)) 1 "^. (2-36) 

We conclude by observing that \x a < C(f3 a ) and using the following standard 
approximation properties of the elements tpi (see for instance |32j). 

jnf II Vu - Vto|U"(n) < Cj(T h )h\\u\\ w 2,2 {n y (2.37) 
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2.3 Galerkin with localized elements. 



Let us prove theorem ll.6l We assume that the coarse mesh is not unadapted 
to F. Let K be a triangle of Th and let a be a node of K such that Tj^ n {K) = 
-Xq where 8 is the interior angle between (F(a), F(b)) and (F(a), F(c)); 
(b, c) being the other nodes of K. Let us prove the following lemma 



,| j— 1 1 1 ct tad 



Lemma 2.3. 

\V F u(K) - V F «(o)| < 3^(A-)||V»| 
Proof. It is easy to check that 

u(6) - u( ) = - F(a))V« o F(a) + (F(b) - F(a)) 

Where the vector q ba is defined by 



Vu[F(a) + s{F{b) - F(a))] - Vu[F(a)] 



ds. 



(2.38) 
(2.39) 

(2.40) 



We will use the notation f ba := (F(b) - F(a))/\F(b) - F(a)\. We will write 
/i~ the unit vector obtained by a 90° rotation of f ba towards f ca . Defining 
q ca as in (|2.4()|) we obtain that 



with 
with 



V F u(K) = S7 F u(a) + k 
k = qba~ ^ft 

fca-i.Qba Qca) 



A 



Which leads us to 

\V F u(K) - V F u(a)\ < 



fca-fba A 

3 



fca ■ fi 



\Vu\ 



c° 



\F\\« a ,h aa 



(2.41) 
(2.42) 
(2.43) 

(2.44) 
□ 



The following lemma is a direct consequence of lemma [^3*1 
Lemma 2.4. Let K € Th and let x G $7 £foen 

|V F n(K) - V F u(x)| < Vfi||c«*(l + \\F\\° a ,)(h + dist(x,K)) ( 

Let us write Z^u the interpolation of u over the space Z^: 



Z h u(x) = ^2 u(xi)£i(x). 



(2.45) 
(2.46) 
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Lemma 2.5. We have 

a*[u-Z h u]<Cr)^\\V^c4 F \\c<*' haa ' D - ( 2 - 47 ) 

Proof. We have 

a K [u-Z h u] = / t (V F u(x)-V F u(K))a(x)(V F u(x)-V F u(K)) dx (2.48) 
Jk 

with o~(x) := *V.FaVi ?1 . Using the change of variables -F(cc) = y we obtain 
that 

ojc[«-2 fc «]= / t (Vu(y)-V F u(K))Q(y)(Vu(y)-V F u(K))dy (2.49) 
from which we deduce that 

a K [u-Z h u] < (3r& in ||V«|M|F||" '^T / sup *eQe. (2.50) 

Jf(E-) |e|=l 

Thus 

a>-£ h u] < C« in || Vtt|| C « 11^1!^^')^ (2.51) 

where -D has been defined by (|2.31l) . □ 

Theorem 11.61 is implied by lemma 12.51 theorem 12.21 lemma 12.21 and the 

following inequality 

a* [u - u f ] < a* [u - Z h u] . (2.52) 
Let us now prove theorem 11.71 By the triangle inequality 

a[u- J h u f ] <a[u- J h u] + a [j h u - J h u f ] . (2.53) 

We write J^u := {Jhu) o F . j^u is a linear interpolation of u on the 
tessellation T F . Now using the identity 

a[u- J h u] = Q[u- J h u] (2.54) 

we obtain that 

a[u - J h u] < \\u\\ c <*h a D. (2.55) 

Where we have written h the maximal length of the edges of T F . Observe 
that h < h a \\F\\ Ca ,. Theorem O 

is a consequence of inequalities (|2.53|) . 
(|2.55|) . lemmas I2.HI and |2,5| theorem 12.21 and the following inequality 

a* [Z^u — u J ~\ < 2a* \u — Z^u^ . 
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Lemma 2.6. 



i 



iJ h u-J h uf] <2S^ir,l^{^^)^a*[Z h u-u^ (2.56) 
and 

a [J h u - J h u f ] < fi a v*a* [Z h u - u f ] (2.57) 
Proof. Let us write w := J^u — Jh u ■ We need to bound a[w]. We have 

a[ w ] = ^^\a*[Z h w]. (2.58) 
a* [Z h w\ 

a*[Zhw] = a*[ZhU — u*\ has already been estimated in lemma l2~51 Observe 
that 

a*[J h w] = Q[w]. (2.59) 
w is piecewise linear on T F . Using property (jl.43j) we obtain that 

Q[i5]<CUQ)^H. ( 2 -60) 

Moreover, observing that 

QoF = r (det(o))' (2.61) 

(det(a)) 2 

we obtain that 

Amax(Q) <4(det(a))5. (2.62) 

Equation (|2.62j) is valid in dimension 2, in dimension higher than 2 we would 
use the inequality (|2.36|) 

1 n-2 1 

A m ax(Q) </4(A min (<x)) 2 (det(a)) 5 . (2.63) 

We now need to bound from below a*[Zhw]/£h[w]. For K & 7^ let us write 
H(K) the matrix 

H(K):= [ t (VF(K))- lt VF(x)aVF(x)(WF(K))- l dx. (2.64) 

We need to estimate infm =1 t lH(K)l. Let us write a, 6, c the nodes of K and 

fix) := - F(a)) {VF{K)) ~ l ■ I. (2.65) 

Let us observe that 

HH{K)l = [ *V/oV/. (2.66) 
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Moreover /(a) = 0, f(b) = (b — a). I and f(c) = (c — a). I. Let us assume 
without loss of generality that \f(b)\/\b — a\ > \f(c)\/\c — a\. Then 

HH(K)l > inf f tVwaVw. (2.67) 

w£C°°(n),w(a)=0,w(b)=f(b) J K 

The quantity appearing in (j2.67j) is the resistance metric distance between 
the points a and b and it is easy to check that (see for instance lemma 1.1 
of 0) 

inf f 'VwoVw > A min (a)Vol(K)(^^) 2 . (2.68) 

w£C™{Q),w(a)=0,w(b)=f(b) J K x \b-a\ / 

Thus 

HH(K)l > A min (a)Vol(K)( l( ^~ Q a )J| ) 2 . (2.69) 
Let us observe that 

\(b-a).l\ 1 
\b-a\ 4r? max 

It follows that 

a*[Z h w] > \\Vl h w\\ 2 L a\ m i n (a) — \ — • (2.71) 



Thus 

^ > W)=3- (2.72) 
which leads us to equation (|2.56|) . 

Remark 2.5. One of the methods employed with non-conformal elements to 
ensure the stability and convergence of the scheme is the so called patch test. 
In our proof the stability condition and convergence is ensured by (|2.72|) and 
an uniform lower bound on ry max . 

To obtain Q2.57|) let us observe that 

a[w]=Q[w\. (2.73) 

Thus 

a[ w ] = V / t Vw{K F )Q{y)Vw{K F )dy. (2.74) 

It follows that 



a[w] <v*J2 f Y^'VwiK^Q^VwiK^dy. (2.75) 
from equation (|2.61() we obtain that 

Amax(Q) 



Amin(Q) 
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< n„. (2.76) 



Next, observing that 

V f t Vw{K F )Q{y)Vw{K F )dy = a*[Z h w]. (2.77) 

T 77n- JF(K) 



i<er h 



we obtain (|2~57|) . □ 

2.4 Numerical homogenization from the informa- 
tion point of view 

In this subsection we will prove theorem 11.81 and theorem 11.91 The method 
introduced in subsection 11.41 can be formulated in the following way: look 
for u m <E Vh such that for all i G Mh, 

a*[<p i ,Z h u m ] = (<p i ,g) I ? a (2.78) 

which implies the following finite volume orthogonality property for all i £ 

a*[ifi,Z h u m -u] = 0. (2.79) 
Let us write w = u — Z^u" 1 . By equation (|1.42|) we obtain that 

/ „ r n \ - 1 a*\v, Zhw] 

(W) 2 <^sup L ' * J . (2.80) 
b «ez h [£ h [v]) 2 

By the orthogonality property we have 

a*[v,Zhw] = a*[v,ZhU — u]. (2-81) 

Thus 

a*[v,Z h w] < (\ max (a))z\\Vv\\ L 2 {n) (a*[Z h u-u\)*. (2.82) 
Using the inequality 

||Vf ||| 2(n) < VmaxShiv]. (2.83) 

We deduce that 

Sh[w]* < -^{KMVm^)Ha*[Z h u - u})K (2.84) 
It follows from ([PH]) that 

\\Vl h u-Vu m \\ L2{n) < ^(X max (a))^(a*[Z h u-u])K (2.85) 
And we deduce from Poincare inequality that 

\\l h u-u m \\ L 2 (n} <C u ^(X max (a))^[a*[2 h u-u])K (2.86) 
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We obtain theorem 11.81 by from equations (|2.85|) . ([2.86ft . lemma 12.51 and 
theorem 12.21 Let us now prove theorem 11.91 Using triangle inequality we 
obtain 

a[u- J h u m ] <a[u- J h u] + a [j h u - J h u m ] . (2.87) 

The object a[u — JhU\ has already been bounded from above by (|2.55|) . 
Writing w := Jh,u — JhU m we have 

= ^fs h {w}. (2-88) 

But £h[w] has already been estimated in equation (|2.84|) . It remains to 
notice that 

a[J h w] = Q[w\. (2.89) 

From this point the arguments are similar to the ones employed in subsection 
E3 indeed 

Q[w] < A max (Q)||W;|| L 2 {n) < v * max X max (Q)£ h [w}. (2.90) 

2.5 Numerical homogenization from a transport 
point of view. 

We assume the mesh to be regular in the following sense: the nodes of the 
Vornoi diagram of 7^ belong elements of the primal mesh. In dimension 2 
this means that for each triangle K £ 7^ the intersection of the median of K 
(the circumcenter) belongs to the interior of K. Let us write Yh the vector 
space spanned by the functions Xi- For v E we define y^v by 

y h v := *iXi- (2-91) 

The metric numerical homogenization method can be formulated in the 
following way: look for u v £ (the space spanned by the elements such 
that for all i G Mh, 

a*[Xi,u v ] = (xi,9)L*n (2-92) 
which implies the following finite volume orthogonality property for all i E 

a*\xi,u v -u]=0. (2.93) 
Equation ()2.92|) can be written 

VV / n.a.V^= f g. (2.94) 
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Let us write w := Z^u — u v . From the equation (|1.58j) we obtain that 

/„ r 1 1 a*\v,w] , 

(S h [w])2< — snp-^-i. 2.95) 

$ v vey h (£ h [v})-2 

Using the orthogonality property of the finite volume method we obtain that 
for v £ y h 

a*[v,w] = - ^2 v i / n.a(X7Z h u- Vu). (2.96) 

Writing £^ the edges of the dual tessellation (edges of the control volumes) , 
we obtain that 



a v, w 



\=^(v j -v i ) n %j .a{VZ h u-Vu). (2.97) 



Where e+j is the edge separating the control volume Vi from the control 
volume Vj and riij is the unit vector orthogonal to eij pointing outside of 
Vi. It follows that 



a*[v,w] <{£h[v\) 2 \\V F Z h u - Vf«||l°°(£*) 
|ejj| 2 A max (a)A max (cr) 



(2.98) 



Now let us observe that 



l%| 2 <3wVol(fi). (2.99) 



(2.100) 



We then have from equation (|2.95|) 

i 3 

i 

Vol(tt)(A max (a)A 

max 

(a)) 2 . 

Equation (|1.59|) of theorem 11.101 is then a straightforward consequence of 
equation (|1.43|) and lemma l2~H Let us now prove equation (|1.60|) of theorem 
11.101 By triangle inequality 

a[u-J h u v ] <a[u- J h u]+a[j h u- J h u v ]. (2.101) 

a[u — Jhii\ has already been estimated in equation H2.55j) . Writing w := 
Jku — Jhu" we have 

= °L^ Eh[w] . ( 2 .i02) 

But £h[w] has already been estimated in equation l2.1()()l It remains to notice 
that g*T^ has already been estimated in equations (|2,59|) and (|2.60j) to 
conclude the proof. 
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3 Numerical Experiments 



Let us now give illustrations of the implementation of this method. The 
domain is the unit disk in dimension two. Equation (jl . lj) is solved on a fine 
tessellation characterized by 66049 nodes and 131072 triangles. The coarse 
tessellation has 289 nodes and 512 triangles (figure EJ). It is important to 
recall that since our methods involve the computation of global harmonic 
coordinates, the memory requirement and CPU time are not improved if one 
needs to solve 1)1- ljl only one time, whereas localized methods such as the 
one of Hou and Wu or E. and Engquist do improve the memory requirement 
or the CPU time. 
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Figure 9: Coarse Grid 

The elliptic operator associated to equation (jl.lj) has been up-scaled to 
an operator defined on the coarse mesh (compression by a factor of 300) 
using 5 different methods: 

• The Galerkin scheme described in subsection 11.21 using the multiscale 
finite element ipi noted FEM_?/>. 

• The Galerkin scheme described in subsection 11.31 using the localized 
elements £j noted FEM_£. 

• The metric based compression scheme described in subsection ll ,4l noted 
MBFEM. 

• The finite volume method described in subsection 11.51 noted FVM. 

• A multi-scale finite element method noted LFEM where F is computed 
locally 5 on each triangle K of the coarse mesh as the solution of a cell 
problem with boundary condition F{x) = x on dK. This method has 
been implemented in order to understand the effect of the removal of 
global information in the structure of the metric induced by F. 

5 instead of globally 
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(a) a. 



(b) r F 



Figure 10: Example ^ Trigonometric multi-scale. 



Example 1. Trigonometric multiscale 

r \ 1 / l.l+sin(27nr/ei) . I.l+sin(27rj//e2) j_ l.l+cos(27nr/e3) . I.l+sin(27rj//e4) . 

a \ x ) — 6U.l+sin(27rj//ei) + l.l+cos(27nr/e 2 ) l.l+sin(27n//e 3 ) + l.l+cos(27nr/<E 4 ) + 
l.l+cos(27ra;/e5) . • (a_J2 2\ i i\ i. 1 1 1 1 

i.i+sin^M) + sm(4xV) + 1), where e x = z ,e 2 = ^, e 3 = ^ ,e 4 = ^ = 

j_ 

65 ' 

Figure ^1 is an illustration of T the deformation of the coarse mesh 
(figure EJ) under the metric induced by F. The deformation is small since 
the medium is quasi-periodic. The weak aspect ratio for triangles the coarse 
mesh in the metric induced by F is r/j^ in = 1.1252. Table^gives the relative 
error estimated on the nodes of the coarse mesh between the solution u of 
the initial PDE (|l.ljl and an approximation obtained from the up-scaled 
operator on the nodes of the coarse mesh. Table |2] gives the relative error 
estimated on the nodes of the fine mesh between u and the ^-interpolation 
of the previous approximations with respect to F on a fine resolution. Figure 



11(a) gives the condition number of the stiffness matrix associated to the 



up-scaled operator versus — log 2 h (logarithm of the resolution). Figure 



11(b) gives the relative Li-distance between u and its approximation on 



the coarse mesh in log scale versus — log 2 h (logarithm of the resolution) . 
Observe that for the method LFEM this error increase with the resolution 
this is an effect of the so called cell resonance observed in |5U| and [2]. This 
cell resonance does not occur with the methods proposed in this paper. The 
finite volume method is characterized by the the best stability and one of 
the best accuracy at a coarse resolution. The increase in the error observed 
for this method as the resolution is decreased is a numerical artifact created 
by the fine mesh: one has to divide the coarse tessellation into coarse control 
volumes. These coarse control volumes are unions of the control volumes 
defined on a fine mesh and when the refinement between the coarse and 
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(a) Condition Number. (b) Coarse Mesh L 1 er- 

ror. 

Figure 11: Example ^ Trigonometric Multiscale 



Coarse Mesh Error 


FEM_^ 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0042 


0.0022 


0.0075 


0.0032 


0.0411 


L 2 


0.0039 


0.0024 


0.0074 


0.0040 


0.0441 


L°° 


0.0059 


0.0090 


0.0154 


0.0117 


0.0496 


H 1 


0.0060 


0.0262 


0.0568 


0.0203 


0.0763 



Table 1: Example ^ Trigonometric multi-scale. 



Fine mesh Error 


FEM_V> 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0042 


0.0085 


0.0053 


0.0080 


0.0593 


L 2 


0.0043 


0.0082 


0.0061 


0.0078 


0.0591 


L°° 


0.0063 


0.0112 


0.0154 


0.0141 


0.0597 


H 1 


0.0581 


0.0540 


0.0778 


0.0601 


0.0943 



Table 2: Example ^ Trigononmetric Multiscale 



the fine mesh is small and the triangulation irregular it is not possible to 
divide the coarse tessellation into control volumes intersecting the edges of 
the primal mesh close to the midpoints of those edges and the other control 
volumes close to the barycenters of the coarse triangles. 
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(a) a. 



(b) r F 



Figure 12: Example 121 High Conductivity Channel. 



Example 2. High Conductivity channel 

In this example a is random and characterized by a fine and long ranged 
high conductivity channel. We choose a(x) = 100, if x is in the channel, and 
a(x) = 0(1), if x is not in the channel. The weak aspect ratio for triangles 
the coarse mesh in the metric induced by F is 7y^ in = 2.2630. Table gives 
the relative error estimated on the nodes of the coarse mesh between the 
solution u of the initial PDE (|1.1|) and an approximation obtained from 
the up-scaled operator on the nodes of the coarse mesh. Table 0] gives 
the relative error estimated on the nodes of the fine mesh between u and 
the ^-interpolation of the previous approximations with respect to F on 
a fine resolution. Figure |13(a)| gives the condition number of the stiffness 
matrix associated to the up-scaled operator versus — log 2 h (logarithm of 
the resolution). Figure [l3(b)| gives the relative Li-distance between u and 
its approximation on the coarse mesh in log scale versus — log 2 h. 

Observe in figure El that the effect of the new metric on the mesh is to 
bring close together nodes linked by a path of low electrical resistance. 

Remark 3.1. Let us recall that the natural distance associated to the Laplace 
operator on a fractal space is also the so called resistance metric |53j . |71| . 
[S]. It is thus natural to find that a similar (not equivalent) notion of dis- 
tance allows the numerical homogenization PDEs with arbitrary coefficients. 
More precisely the analogue of the resistance metric here are the harmonic 
mappings. The analysis of these mappings allows to bypass boundary layers 
effects in homogenization in periodic media it allows to obtain quanti- 
tative estimates on the heat kernel of periodic operators (i6 or to analyze 
PDEs characterized by an infinite number of non separated scales |14j . |67| . 
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(a) Condition Number. (b) Coarse Mesh L 1 er- 

ror. 

Figure 13: Example High Conductivity Channel 



Coarse Mesh Error 


FEM_^ 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0022 


0.0081 


0.0127 


0.0062 


0.0519 


L 2 


0.0025 


0.0096 


0.0179 


0.0081 


0.0606 


L°° 


0.0120 


0.0227 


0.0549 


0.0174 


0.1223 


H 1 


0.0120 


0.0384 


0.0919 


0.0265 


0.1514 



Table 3: Example |2j High Conductivity Channel. 



Fine mesh Error 


FEM_ip 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0070 


0.0155 


0.0164 


0.0121 


0.0612 


L 2 


0.0069 


0.0153 


0.0202 


0.0123 


0.0743 


L°° 


0.0133 


0.0227 


0.0573 


0.0214 


0.1226 


H 1 


0.0760 


0.1032 


0.1838 


0.0820 


0.2142 



Table 4: Example El High Conductivity Channel. 
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(a) a. 



(b) r F 



Figure 14: Example El Random Fourier Modes . 



Example 3. Random Fourier modes. 
In this case, a(x) = e h ( x \ with 

h(x) = ^2 ( a k sin(27r/c • x) + bk cos(27rfe • x)) 
\k\<R 

where and bk are independent identically distributed random variables 
on [—0.3, 0.3] and R = 6. This is an other example where scales are not 
separated. The weak aspect ratio of the triangles in the metric induced by 
F is rfL^ = 3.4997. Observe the deformation induced by the new metric 
(figure fTl]) . Observe that distances between u and the interpolation of the 
coarse mesh approximations to the fine mesh are larger (tables El and EJ), this 
is due to the fact that those errors depend on the aspect ratio ^ (which 
is not the case for the coarse mesh errors) of course one could improve the 
compression by adapting the mesh to the new metric but this has not been 
our point of view here. We have preferred to show raw data obtained with 
a given coarse mesh. The figures [TBI and IT71 give the L 1 , L 2 , L°° and H 1 
relative error (log 2 basis versus log 2 basis of the resolution). The x-axis 
corresponds to the refinement of coarse mesh, the y-axis is the error. The 
tables □ and H give the convergence rate in different norms (the parameter 
a in the error of the order of h a ). 
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Coarse Mesh Error 


FEM_^ 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0027 


0.0075 


0.0117 


0.0106 


0.1197 


L 2 


0.0028 


0.0087 


0.0130 


0.0125 


0.1169 


L°° 


0.0066 


0.0278 


0.0320 


0.0376 


0.1358 


H 1 


0.0133 


0.0648 


0.0805 


0.0597 


0.1514 



Table 5: Example El Random Fourier Modes 



Fine mesh Error 


FEMj0 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0112 


0.0148 


0.0148 


0.0188 


0.1304 


L 2 


0.0177 


0.0223 


0.0184 


0.0202 


0.1265 


L°° 


0.0773 


0.0824 


0.0614 


0.0680 


0.1669 


H 1 


0.0972 


0.1152 


0.1307 


0.1659 


0.1725 



Table 6: Example El Random Fourier Modes. 



Method 


L 1 


L 2 


L°° 


H 1 


FEM_V> 


1.62 


1.66 


1.56 


1.44 


FEM_£ 


1.38 


1.27 


1.23 


1.18 


MBFEM 


1.38 


1.40 


1.27 


1.08 


FVM 


0.53 


1.14 


1.26 


1.03 


LFEM 


1.51 


1.53 


1.62 


1.46 



Table 7: Coarse mesh approximation convergence rate 



Method 


L 1 


L 2 


L°° 


H 1 


FEM_ip 


1.74 


1.61 


1.23 


0.89 


FEM_£ 


1.57 


1.47 


1.23 


0.91 


MBFEM 


1.54 


1.52 


1.21 


0.96 


FVM 


0.75 


1.16 


1.22 


0.58 


LFEM 


1.52 


1.54 


1.42 


1.10 



Table 8: Fine mesh approximation convergence rate 
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4.5 5 



4.5 5 



(a) Condition Number. 



(b) Coarse Mesh L 1 er- 
ror. 



Figure 15: Example El Random Fourier Modes 
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(a) L 1 error. 



(b) L 2 error. 




(c) L°° error. (d) H 1 error. 

Figure 16: Coarse mesh error (log 2 ) L 1 , L 2 , L°° and H 1 errors v.s coarse mesh 
refinement, Example EJ Random Fourier Modes. 
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(a) L 1 error. 



(b) L 2 error. 




45 



1 1 -1 -0.5 0.5 

(a) a. (b) T F 

Figure 18: Example EJ Random Fractal. 



Coarse Mesh Error 


FEM_^ 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0024 


0.0075 


0.0231 


0.0073 


0.0519 


L 2 


0.0025 


0.0085 


0.0241 


0.0100 


0.0606 


L°° 


0.0094 


0.0399 


0.0920 


0.0398 


0.1694 


H 1 


0.0161 


0.0718 


0.1553 


0.0493 


0.3107 



Table 9: Example IH Random Fractal. 



Example 4. Random fractal 

In this case, a is given by a product of discontinuous functions oscillating 
randomly at different scales, a{x) = a\{x)a>2,(x) ■ ■ ■ a n (x), and ai{x) = c pq 
for x G [t|-, ^r-) x [^-, ^r-), c pq is uniformly random in [^,7], n = 5, and 
7 = 2. The weak aspect ratio is r/^ in = 2.4796. Table |H1 gives the relative 
error estimated on the nodes of the coarse mesh between the solution u of 
the initial PDE an d an approximation obtained from the up-scaled 

operator on the nodes of the coarse mesh. Table ITUl gives the relative error 
estimated on the nodes of the fine mesh between u and the ^-interpolation 
of the previous approximations with respect to F on a fine resolution. Figure 



19(a) gives the condition number of the stiffness matrix associated to the up- 
scaled operator versus — log 2 /i (logarithm of the resolution). Figure [l9(b)| 
gives the relative Li-distance between u and its approximation on the coarse 
mesh in log scale versus — log 2 h. 
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(a) Condition Number. (b) Coarse Mesh L 1 er- 

ror. 

Figure 19: Example Random Fractal 



Fine mesh Error 


FEM_?/> 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0108 


0.0147 


0.0245 


0.0142 


0.0765 


L 2 


0.0155 


0.0198 


0.0280 


0.0173 


0.0812 


L°° 


0.0662 


0.0802 


0.0919 


0.0720 


0.1694 


H 1 


0.1015 


0.1231 


0.1838 


0.1433 


0.2642 



Table 10: Example EJ Random Fractal. 
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(a) a. (b) T F . 

Figure 20: Example El Percolation. 



Example 5. Percolation at criticality 

In this case, the conductivity of each site is equal to 7 or I/7 with proba- 
bility 1/2. We have chosen 7 = 4 in this example. Observe that some errors 
are larger for this challenging case because a percolating medium generates 
flat triangles in the new metric, indeed r/^ in = 22.3395. Table fTTI gives the 
relative error estimated on the nodes of the coarse mesh between the solution 
u of the initial PDE and an approximation obtained from the up-scaled 
operator on the nodes of the coarse mesh. Table IT2l gives the relative error 
estimated on the nodes of the fine mesh between u and the ^-interpolation 
of the previous approximations with respect to F on a fine resolution. Fig- 
ure 21(a) gives the condition number of the stiffness matrix associated to 



the up-scaled operator versus — log 2 h (logarithm of the resolution). Figure 



21(b) gives the relative Li-distance between u and its approximation on the 
coarse mesh in log scale versus — log 2 h. Observe that the methods based 
on a global change of metric do converge but when that numerical homog- 
enization is done by computing only local coarse parameters (in averaging 
or finite elements techniques) then convergence is not guaranteed without 
further assumptions on a (see the curve of LFEM). 
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(a) Condition Number. (b) Coarse Mesh L 1 er- 

ror. 

Figure 21: Example Percolation 



Coarse Mesh Error 


FEM_^ 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0034 


0.0253 


0.0485 


0.0167 


0.2848 


L 2 


0.0041 


0.0265 


0.0523 


0.0189 


0.2851 


L°° 


0.0163 


0.0813 


0.0643 


0.0499 


0.3018 


H 1 


0.0343 


0.0843 


0.1070 


0.0713 


0.3740 



Table 11: Example Percolation. 



Fine mesh Error 


FEMjip 


FEM_£ 


MBFEM 


FVM 


LFEM 


L 1 


0.0115 


0.0265 


0.0585 


0.0216 


0.3024 


L 2 


0.0152 


0.0268 


0.0628 


0.0229 


0.3015 


L°° 


0.0500 


0.0527 


0.0940 


0.0497 


0.3135 


H 1 


0.1000 


0.1712 


0.1954 


0.1343 


0.3964 



Table 12: Example Percolation. 
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Coarse Mesh Error 


FEM_ipi in 


FEM_^ sp 


L 1 


0.0437 


0.0046 


L 2 


0.0426 


0.0052 


L°° 


0.0614 


0.0096 


H 1 


0.0746 


0.0227 



Table 13: Example HH Random Fourier Modes 



3.1 Numerical experiments with splines. 

We have seen that if a is stable then u o F _1 belongs to W 2,p (£l) with 
p > 2. It is thus natural to expect a better accuracy by using C 1 -continuous 
elements in the method described in subsection 11.21 instead of piecewise 
linear elements. This increase of accuracy has already been observed in |2j 
when F is obtained as the solution of a local cell problem. In our case 
(the harmonic coordinates are computed globally) we also observe a sharp 
increase of the accuracy of finite element method of subsection 11.21 bv using 
splines for the elements tp^. 

We refer to |48| I24j for methods using C 1 finite element. One possibility 
is the weighted extened B-splines (WEB) method developed by K. Hollig in 
[H1[1B], these elements are C 2 -continuous. They are obtained from tensor 
products of one dimensional elements. The Dirichlet boundary condition is 
satisfied using a smooth weight function uj, such that u = at the bound- 
ary. The condition number of the stiffness matrix is bounded from above 
by 0(h~ 2 ) (we have the same optimal bound on a Galerkin system with 
piecewise linear elements). 

We have considered two challenging multi-scale medium for our numeri- 
cal experiments: random Fourier modes and percolation. For the simplicity 
of the implementation a square domain has been considered, and weighted 
spline basis are used instead of the WEB spline basis. For a square domain 
[—1, 1] x [—1, 1], the weight is u = (1 — x 2 )(l — y 2 ). Two methods have been 
compared, 

• The Galerkin scheme using the finite elements fa = <pi o F, where ipi 
are the piecewise linear nodal basis elements of subsection II. 2[ noted 
FEM_fa in 

• The Galerkin scheme using the finite element fa = (pi o F, where fa 
are weighted cubic B-spline basis elements, noted FEMj0 sp 

The error obtain with the method FEM_0 sp is much smaller than the 
one obtained with the method FEM_^ n as it is shown in figures [2*2*1 and I2U1 
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Fine mesh Error 


FEMj0 lin 


FEM_^ sp 


L 1 


0.0546 


0.0077 


L' 2 


0.0529 


0.0096 


L°° 


0.0920 


0.0289 


H 1 


0.2109 


0.0547 



Table 14: Example H3 Random Fourier Modes. 



Coarse Mesh Error 


FEM_1pii n 


FEM_^ sp 


L 1 


0.0393 


0.0080 


L' 2 


0.0379 


0.0098 


L°° 


0.0622 


0.0309 


H 1 


0.0731 


0.0404 



Table 15: Example Percolation 



Fine mesh Error 


FEM_ipi in 


FEM_^ sp 


L 1 


0.0470 


0.0099 


L' 2 


0.0464 


0.0130 


L°° 


0.1174 


0.0554 


H 1 


0.2030 


0.0838 



Table 16: Example El Percolation. 
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(a) Condition Num- (b) Coarse Mesh L 1 (c) Fine Mesh L 1 er- 

ber. error. ror. 

Figure 22: Example El Random Fourier Modes. 
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(a) Condition Num- (b) Coarse Mesh L 1 

ber. error. 




2 2.5 3 3.5 i 4.5 5 



(c) Fine Mesh L 1 error. 
Figure 23: Example 03 Percolation at criticality. 
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